function [ Price ] = UncertainVol(S,Asset,sigmamax,sigmamin,r,T,N,K1,K2,scenario,theta )
%UNCERTAINVOL under the uncertain volitalitymodel we determine the price of
%the butterfly spread with strikes K1,K2,K3
num_it=0;
q=numel(S);
delta_t=T/N;
M= spalloc(q,q,3*(q-2)+1);
Large=10^6;
tol=1/Large;
I=speye(q);
Vn=max(S'-K1,0)-2*max(S'-(K1+K2)*0.5,0)+max(S'-K2,0);
[delta gamma]=greeks(S,Vn);
sigma=getsigma(sigmamax,sigmamin,gamma,scenario);
M(1,1)=r*delta_t;

% Test
A1 = zeros(q-2,1);
B1 = zeros(q-2,1);
for j=2:q-1
    [alpha beta]=PCFDDcoeff(S,j,r,sigma);
    M(j,j-1:j+1)=delta_t*[-alpha,(alpha+beta+r),-beta];
    A1(j-1) = alpha;
    B1(j-1) = beta;
end
Mk1=M;
V2 = Vn;
for n=1:N
    const=1+tol;
    Vn1=zeros(q,1);
    while (const>tol)
        num_it=num_it+1;
        Vn1=(I+(1-theta)*Mk1)\((I-theta*M)*Vn);
        [delta gamma]=greeks(S,Vn1);
        sigma=getsigma(sigmamax,sigmamin,gamma,scenario);
        for j=2:q-1
            [alpha beta]=PCFDDcoeff(S,j,r,sigma);
            Mk1(j,j-1:j+1)=delta_t*[-alpha,(alpha+beta+r),-beta];
        end
        const=max(abs(Vn1-V2)./max(1,abs(Vn1)));
        V2=Vn1;
    end
    Vn=Vn1;
    M=Mk1;
end
Price=Vn(find(S==Asset));
greekplot(S,gamma,delta,Vn,theta);
end


function[alpha,beta]=PCFDDcoeff(S,j,r,sigma)
alpha_c=((sigma(j-1)^2)*S(j)^2)/((S(j)-S(j-1))*(S(j+1)-S(j-1)))-r*S(j)/(S(j+1)-S(j-1));
beta_c=((sigma(j-1)^2)*S(j)^2)/((S(j+1)-S(j))*(S(j+1)-S(j-1)))+r*S(j)/(S(j+1)-S(j-1));

alpha_f=((sigma(j-1)^2)*S(j)^2)/((S(j)-S(j-1))*(S(j+1)-S(j-1)));
beta_f=((sigma(j-1)^2)*S(j)^2)/((S(j+1)-S(j))*(S(j+1)-S(j-1)))+r*S(j)/(S(j+1)-S(j));

if(alpha_c>0&&beta_c>0)
    alpha=alpha_c;
    beta=beta_c;
elseif(alpha_f>0&&beta_f>0)
    alpha=alpha_f;
    beta=beta_f;
else
    alpha=((sigma(j-1)^2)*S(j)^2)/((S(j)-S(j-1))*(S(j+1)-S(j-1)))-r*S(j)/(S(j)-S(j-1));
    beta=((sigma(j-1)^2)*S(j)^2)/((S(j+1)-S(j))*(S(j+1)-S(j-1)));
end

end
function[delta gamma]=greeks(S,V)
S=S';
deltaS1=S(3:end)-S(2:end-1);
deltaS_1=S(2:end-1)-S(1:end-2);
V1=V(3:end)-V(2:end-1);
V2=V(2:end-1)-V(1:end-2);
V3=V(3:end)-V(1:end-2);
delta=V3./(deltaS1+deltaS_1);
gamma=((V1./deltaS1)-(V2./deltaS_1))./(0.5*(deltaS1+deltaS_1));
end
function[sigma]=getsigma(sigmamax,sigmamin,gamma,scenario)
if scenario=='w'
    sigma=sigmamax*(gamma<=0)+sigmamin*(gamma>0);
else
    sigma=sigmamin*(gamma<=0)+sigmamax*(gamma>0);
end
end
function greekplot(S,gamma,delta,V,theta)
ind50=find(S==50);
ind150=find(S==150);
subplot(1,3,1);
plot(S(ind50:ind150),V(ind50-1:ind150-1));
if (theta==0)
    title('Price of European Put as a function of Stock price under CN-Ranacher Method with constant TimeSteps');
else
    title('Price of European Put as a function of Stock price under CN-Ranacher Method with variable TimeSteps');
end
xlabel('Stock Price');
ylabel('Option Price');
subplot(1,3,2);
plot(S(ind50:ind150),delta(ind50-1:ind150-1));
if (theta==0)
    title('Delta of American Put as a function of Stock price under CN-Ranacher Method with constant TimeSteps');
else
    title('Delta of American Put as a function of Stock price under CN-Ranacher Method with variable TimeSteps');
end
xlabel('Stock Price');
ylabel('put delta')
subplot(1,3,3)
plot(S(ind50:ind150),gamma(ind50-1:ind150-1));
if (theta==0)
    title('Gamma of American Put as a function of Stock price under CN-Ranacher Method with constant TimeSteps');
else
    title('Gamma of American Put as a function of Stock price under CN-Ranacher Method with variable TimeSteps');
end
xlabel('Stock Price');
ylabel('put gamma');

end
